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i. INTRODUCTION 


Spectral measurements are frequently required in fluid mechanics 
applications. Traditionally they have been made using analog techniques. 
With the development of the Fast Fourier Transform algorithms in the 
mid 1960's, digital techniques have evolved which enable power spectral 
densities and correlation functions to be calculated with costs much 
less than were previously possible. This report is intended to describe 
tive Paces: Terrie, Taeintyas eigouluuas avallenis of Ceterpte Seas 
University, outline some of the difficulties encountered in using these 
algorithms, and provide a brief description of actual computer programs 


being used for spectral analysis on the CDC 6400 computer. 
2. EXISTING FFT ROUTINES AT CSU 


There are presently a number of computer programs used at CSU 
which use available FFT routines. Two FFT programs are being used 
extensively by the Fluid Mechanics and Wind Engineering group. These 
are FOR2D and FOURT, both a part of the IBM Contributed Program Library. 
FOURT (IBM Contributed Program No. 360D-13.4001) is presently on the 
system Fortran library (FTNLIB). FOR2D (IBM Contributed Program 
No. 360D-13.4006) is usually stored on a permanent file. For CSU users, 
a deck or access to this permanent file may be obtained by contacting 
Robert Akins or Dr. J. Peterka. The major difference between these two 
programs is that FOURT is written to use data located in the core of the 
computer and FOR2D is written to use data located on an external storage 
device. More detailed comments on these two specific subroutines appear 


in later sections. 


Si: SINGLE CHANNEL FORWARD/ INVERSE TRANSFORM 


Two separate uses of the FFT will be described; (1) calculation 
of a power spectral density from a time series and (2) transformation 
of a power spectral density to obtain an autocorrelation function. An 
explanation of the details of these types of caicniarions can be found 
in Bendat and Piersol (1). It will be assumed in the following dis- 
cussions that the reader is familiar with this reference or an 
equivalent text, 

The single-channel forward/inverse transform is perhaps the most 
Straightforward application of the FFT, and is a good starting point 
for someone beginning to work with the FFT. A useful exercise is to 
select a known fourier transform pair and to perform the same transform 
using the FFT, An example utilizing this type approach will be dis- 
cussed in order to illustrate usage of subroutine FOURT. Appendix Al 
contains a program listing of subroutine,FOURT. A short section of 
comments appears at the beginning of the listing and explains the calling 
parameters and some basic aspects of usage. Use of the program can be 
understood without a detailed understanding of the details of the 
program itself. 

The example transform pair to be used consists of R(t) = a, 

t > 0 and its inverse fourier transform G(w) = 4/1+w*, w > 0. Such 
an R function is often used to represent the autocorrelation function 
of a fluctuating velocity signal and is not only an easy function to 


deal with, but also is of some physical significance. A sample. 


program (Program CHECK) which was written to take a forward and inverse 


~~ a 


LE — 


fourier transform is listed in Appendix Bl. The following discussion 


will be based upon output from that program. References to the 
program will be by line number of the listing in the appendix. 


The program was written to calculate a selected number of values 


A 


This array of values of (R(t), called (D')in the program, is reflected 


prior to performing the forward transform. This reflection is an. 


important operation which is not adequately discussed in most texts. | 


It is needed to satisfy continuous, even-function characteristics of 


of the function R(t) at a time step specified by an input parameter. 


the transform. In using a digital transform technique, one assumes 


that the data record is of infinite length. In order to create a 


record which resembles an infinite record, the function to be trans- 
——— Ai hit 
formed _is reflected about its endpoint, creating a symmetric, even, 
continuous function. A schematic of this reflection and the resulting 
periodic function is shown in (Figure 1. (note that the function is one 
Ne ; ee re gs 
AT increment short of returning to the zero-time value of one at the 
time position 2N AT. This occurs because the reflection point is 


really at the N(AT+1) point rather than precisely on the NAT point 


aS might be expected. The effect of not reflecting R prior to the 


—. 


‘™ 
transform is shown i Figure 2. In other words, the reflection should 


be done about a zero lag time, such that R(t) = R(-tT) so that the 


ee 


reflected correlation function is even. If there are 2N total 


points, N prior to reflection, then R(N+I) = R(N+2-1) for 2 < I < N. 


This scheme of reflection will result in the point R(N+l) not being 
defined._-Since a correlation is normally small at the maximum lag 


it is easiest tQ 1 N+1) equal R(N)\ The reflection in 


program check is performed in lines 35-40. 


_ — — — — — 


After the data has been reflected and the D(2,1)'s set to zero, 
subroutine FOURT is called in line 43. It is necessary to create D- 
as a two-dimensional array because the input to FOURT is complex. In 
many applications of these techniques only real signals are dealt with 
and the complex arithmetic capabilities of the program are not used. 

In these cases the fifth calling parameter of FOURT is set equal to 
zero indicating a real input only. It is important to understand the 
difference between NUMBER and NUMBE2. NUMBER as used in the program 
is the number of points in the unreflected R_ function. NUMBE2 is 
_twice NUMBER or the length of the reflected R _ function. 
After FOURT has been called, the output must be multiplied by a 
_scaling factor in order to obtain the correct G. This factor for 
FOURT is 2 * DELTAT where DELTAT is the timestep of R. This multiplica- 
tion is carried out in line 52 of Program Check. The G_ function 
returned from FOURT has data. uniformly spaced in frequency with the 
data points at | f = n/ (NUMBER * DELTAT),] n = 0, NUMBER-1. 


Prior to performing the inverse transform, G is also reflected. 


In addition to reflecting G, the imaginary part of the D array is 
set to zero. Normally small values will appear in this position of 
the array during a transform and the inverse transform will be more 
waocastts if the imaginary (D(2,1)) part of the D array is forced to 


zero. These operations are carried out in lines 57 to 65. 


An inverse transform is performed to obtain the original R(t). 
oe ee ee 


Again the output of FOURT must be multiplied by a constant. For the 


inverse transform, this factor is| 1/(4*NUMBER*DELTAT)] The product 


of the factors for the forward and inverse transforms is 1/(2*NUMBER) 


or 1/NUMBER2. This agrees with the factor given in the write-up of 


——<—<————— 
= 


FOURT in Appendix Al, — 


Several examples using the test functions R and G_ will now be 


discussed. A number of different experiments were conducted to determine 


the effect of the total time and the time increment on the accuracy of 


the results. In the tables and plots that forlow (RA will denote the 


_recovered R_ function after a forward and an inverse transform. 


Figure 3 lis a comparison of G for various values of NUMBER, the 

total number of points, and DELTAT, the time step between points. 
These variables were selected to result in all of the R's _ being 
defined for the same total time. This plot is only for the higher 
frequencies; below w = 3, all of the cases tested agree. The N's 
on the plot are for the unreflected data. It can be seen that as 
NUMBER increases and DELTAT decreases the region over which the 
transform is accurate increases. For NUMBER equal to 2048, the results 
are very close to the actual function G for the entire range plotted. 
At higher frequencies, even the case for NUMBER equal to 2048 will 
deviate from the actual values. 

A comparison of _R* for the same cases is shown in Figure 4. 


R* is the initial function R after having been subjected to both 


a forward and inverse transform. In this case virtually all values 


of NUMBER yield an acceptable value for R%*. 
The central processor (Cp) times required for the various 


values of NUMBER2 are shown in Figure 5. It can be seen that there is 


an almost linear increase in Cp time with increasing NUMBER2. These 
f 


times are for the actual transform only; any multiplication or other 


manipulations with the data would increase them. These test cases 
were run under the Scope 3.3.14 on the CSU CDC 6400 computer. 

In summary, the FFT can be used rapidly and economically to 
perform a digital fourier transform of known data. fare must be taken 


ie 


to insure the data are reflected properly prior to the transform and 
i Se a ee = a 


that the appropriate factors are used after the transform. A user with 


—— aa — — 
a a — 


no experience with FFT is strongly urged to experiment with this type 
of application prior to attempting to obtain a spectrum directly from 


digital data. 
4, CALCULATION OF A POWER SPECTRAL DENSITY FROM A TIME SERIES 


Another valuable application of the FFT is the calculation of a 
power spectral density function from a time series. A detailed 


explanation of this process is given in Chapter 9 of Bendat and 


Piersol (1) or in Chapter 6 of Enochson and Ontes (2). The basic 
— ~— —=4 edi i nani. ; 


equations used are straightforward and apply to any FFT routine. 


here is gone significant difference between the procedures outlined 
in these references and the procedure recommended in this report. 


This difference has to do with the addition of zeros to the 


initial time series to avoid having a distorted autocorrelation 
i =e 


at 


function as discussed on pages 312-314 of Bendat and Piersol (1). if 
a = SS ee SSS Sa ad ——. 


one uses the reflection techniques described in Section 3 in obtaining 


an autocorrelation function from a power spectral density, the 


addition of zeros to the initial time series is_unnecessary> This 
ee _ rs $$$ __ 
results in a significant advantage in that the same time series can 
be placed in a data array half the size required if the technique 


described in Bendat and Piersol (1) is used. In other words if a time 


series of data consisting of 2000 points was to be examined using 
Standard procedures, an array of length 4000 would be required. If 
the reflection technique was used, the required array length would 
only be 2000 and there would be a savings in core of 50 percent. In 
most cases the size of the data array is limited by the available core 
of the computer, and therefore the ability to use a smaller array can 
be a significant advantage. 

At this point nomenclature comparable to that in Bendat and 
Piersol will be introduced to make the following discussions easier 
"£6 follow. Denote the time series by x(t), n= 1,N, the fourier 
transform of this time series by X(f, N) n = 1, N, and the spectral 
density function of the time series x by cme aa n=1,N. fT = 
(n-1)/T where T =N At. At is the time increment of the initial time 
series. 


——— 


In terms of these variables, a technique for computation of 


power spectral densities is: 


1. Truncate the data sequence or_add zeros such that N isa 
power of 2. In most cases, the data should be taken to 
provide N data values without adding any zeros. 


2. Yaper this sequence using a cosine taper window. This process 
is discussed in Bendat and Piersol (1), pp. 322-324. 


3. Compute X(t, N) using a FFT routine. 
4. Compute G (f_) using the equation 6 (£) een k. x : 
xn > aa «aes Sa N | K | 


5. Smooth G (£) using either frequency Or segment averaging. 


Frequency averaging averages together severa] values of G, 


from one transform about some value . and replaces all values 
averaged with one average value. Segment averaging is an ensemble aver- 
——— 


age at each sane OF oy of a number of separate transforms. 


These steps are the basis for two programs which will be used as 
examples. It should be noted that a real data sequence Xn will have 
a complex fourier transform X(f N). In a sense the real part is the 
coefficient of the cosine term and the imaginary part is the coefficient 
of the sine term. Therefore in step (4) when the power spectral density 
estimate is computed, the sum of the square of these two terms is used. 
normalized with the variance of the time series. This normalized 
power spectral density will be called PCa or F(n). 

The smoothing in step (5) 1s one of the more subjective aspects 
of the procedure and the technique used will depend upon the type of 
Signal being analyzed, the amount of computer time available, and the 
final use of the power spectral density. The smoothing and the choice 


of N and AT will determine the frequency range of the smoothed 


————, 


nes ee ee There is some choice available in the deter- 
neo of these parameters, and this choice should be made prior to 
taking the data. 

The largest value of N_ which can be used in core with the 
CDC 6400 is 8192 (2°), This is the largest power of two which can 
be used for a data array and not exceed the available core. Frequencies 
will then run from 0 to (F-ek. But T = NAT, and therefore the 
frequencies will run from 0 to F-l*ae . For large N this is 
approximately 1/2AT, the Nyquist frequency. The zero frequency value 
is generally net reliable because the record lengths are of a finite 
length. If the value were to be nearly exact, the total time of the 
input data record, T, should approach infinity. The increment between 


me 


points is equal to - where T is the length in time of the input 


data record. Recall T is equal to NAT. Therefore the high 
frequency end of the power spectral density is determined by AT, the 
time interval of the data record, and the low frequency end is deter- 
mined by the length of the data record. 

Normally the type signal to be examined will dictate the sample 
rate, 1/AT. Once this is determined, and if the maximum range possible 
is desired, N is 8192, the low frequency end of the power spectral 
density is also set. The following table gives these limits for 
sample rates available on the Systems-Development A-D system currently 


in use. 


TABLE 1 - LIMITS FOR POWER SPECTRAL DENSITY COMPUTATION - SYSTEMS - 
DEVELOPMENT A-D SYSTEM. RECORD LENGTH = 8192. 


SAMPLE LOWER UPPER 
RATE LIMIT LIMIT 
AT (SEC) (1/SEC) (HZ) (HZ) 
004 250 031 125.0 
002 500 061 250.0 
001 1000 122 500.0 
0005 2000 244 1000.0 
.00025 4000 .488 2000.0 


If a smaller range is desired, N may be reduced and there will be 
a savings in computer costs. If a larger range is desired, a program 
is available which allows larger N's to be used by employing an 
external storage device such as a disc. This program will be discussed 


later in this section. 
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Smoothing of the power spectral density is required. Two 
techniques are available: segment averaging and frequency averaging. 
These may be used independently or in a combined manner. In segment 
averaging, a number of power spectral densities are computed from 
separate records from the same signal. These estimates of the power 
spectral densities are treated as an ensemble, and an ensemble average 
computed. The number of segments used is determined by the quality of 
the smoothed power spectral density desired and the amount of computer 
time to be expended. Segment averaging will not alter the frequency 
range of the power spectral density, the upper limit will be 1/2AT 
and the lower limit 1/T. 

Frequency averaging involves averaging adjacent points of the power 
Spectral density estimate from one data record. For example every Mm 
points could be averaged and replaced by one point at the midpoint of 
the frequency range of the original m points. This type of averaging 
will have a negligible effect on the high frequency limit of the power 
spectral density, but will normally raise the lower limit substantially, 
depending, of course, on the choice of m and the original AT. 

Factors which enter into the choice of frequency smoothing 
techniques are determined by the ultimate use of the power spectral 
density. If a well-smoothed plot is the desired output, a combination 
of frequency and segment averaging may be employed. If a correlation 
function is to be computed from the power spectral density function, 
then equal frequency spacing must be preserved. Also, the time spacing 
of the correlation obtained is determined by the frequency interval of 
the power spectral density and this relationship should be considered 


in any frequency smoothing. 
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4.1 Calculation of Power Spectral Densities Using Segment Averaging 
Techniques 


In order to provide some examples of the use of both the FFT and 
the averaging techniques, output from a specific program will be 
presented. This program, SEGEMNT, is listed in Appendix B2, and 
references will again be made to line numbers in the program. 

This program follows the suggested routine for computation of 
a power spectral density. Lines 107-133 read one block of data 8192 
elements long off of the data tape, tape 1, and compute the mean and 
the rms of that data record. Lines 138-151 remove the mean from the 
data and divide by the rms to obtain a rms of 1.0. This section of 
the program also tapers the data. Lines 156-167 perform a forward 
fourier transform of the array D, and segment average into array 
SEGMEN. Lines 171-194 reflect the segment averaged spectra and perform 
an inverse transform to obtain a correlation function. The remainder 
of the program is concerned with output and plots of both the correla- 
tion and the power spectral density. Frequency averaging is performed 
in lines 246-260. 

Some sample results from this program will now be used to illustrate 
the effect of segment and frequency averaging. Segment averaging can 
be evaluated using both qualitative and quantitative methods. The 
appearance of both the smoothed spectra and the autocorrelation can 
be compared for different numbers of segments. Figure 6 shows four 
different segment averaged spectra computed from the same data record. 
All four of these spectra were also smoothed using Remi: auccmadinn 
over the high frequency portion. The portion of the title which is of 
the form xx-8192 indicates how many segments of length 8192 were used 
in the calculation of the spectra. It can be easily seen that as the 


total number of records increases, the spectra become smoother. If 
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the spectra are compared by laying one on another, there is no change 

in the best line that could be drawn through the data. In other words, 
if the 64-8192 case is compared with the 4-8192 case, the mean curves 
are identical. Additional qualitative comparisons can be made using 

a number of different criteria. The effective bandwidth, number of 
degrees of freedom and normalized standard error for the different cases 
can be computed using equations (9.140) to (9.149) of Bendat and 

Piersol (1). These values for the cases plotted in Figure 6 are shown 
in Table 2. As the number of segments averaged increases, the 
normalized standard error decreases. The effect of frequency averaging 
in reducing the normalized standard error can also be seen. Another 
means of comparison is available in terms of more physically relevant 
parameters. The area under the spectrum is compared in Table 3 for the 
four cases shown in Figure 6. There is very little difference in these 
integrated quantities as the number of segments increases. These values 
were all computed for the segment averaged spectra before frequency 
averaging. Close attention Should be paid to the integral of F(n). A 


value which is not very close to 1.00 is an indication that, for some 


reason, an incorrect spectrum has been obtained. 

Figure 7 shows the qualitative effect of frequency averaging. All 
three cases were averaged over the same number of segments, and the 
differences are a result of frequency averaging alone. The last line 
of the titles indicate the type of frequency averaging used. The 
different averaging schemes are: (1) no frequency averaging (2) HF 


AVG 10 - no frequency averaging from 0-5.98 HZ, 10 points averaged 


13 


TABLE 2 - EFFECT OF SEGMENT AVERAGING ON POWER SPECTRAL 
DENSITY ESTIMATES 


NUMBER RANGE 
OF POINTS OF EFFECTIVE NUMBER NORMALIZED 
FREQUENCY | SMOOTHING BANDWIDTH OF DEGREES '— STANDARD 
CASE AVERAGED (HZ) (HZ) OF FREEDOM ERROR 
4-8192 1 0-5.98 123 8 500 
10 5.93-500.0 1.220 80 158 
8-8192 1 0-5.98 123 16 353 
10 5.98-500.0 1.220 160 ‘112 
16-8192 1 0-5.98 122 32 250 
10 5.98-500.0 1.220 320 079 
64-8192 l 0-1.09 192 128 125 
3 1.09-5.98 . 366 384 io72 
10 5.98-500.0 1.220 1280 039 


TABLE 3 - COMPARISON OF INTEGRATED PROPERTIES OF POWER 
SPECTRAL DENSITY ESTIMATES 


500 
CASE | F(n)dn 
O 
4-8192 1.014 
8-8192 1.004 
16-8192 1.002 


64-8192 1.007 
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from 5.98-500 HZ (3) HF AVG 10 LF AVG 3 - no frequency averaging from 
0-1.098 HZ, 3 points averaged from 1.098-5.98 HZ and 10 points averaged 
from 5.98-500 HZ. Table 4 is comparable to Table 2 and shows the 

a aia iam 
effective bandwidths, number of degrees of freedom and normalized 


« a pr er ae 
i 


standard error for the cases shown in Figure Ge These criteria are 


ane 


—— 


the only ways to evaluate frequency averaging. In most cases frequency 
averaging will be used to provide a smooth plot of the spectra, and 
the means of frequency smoothing selected will be dependent upon the 
type of data being considered, the frequency range of interest, and the 


ultimate use of the plot. 


TABLE 4 - EFFECT OF FREQUENCY AVERAGING ON POWER SPECTRAL 
DENSITY ESTIMATES 


NUMBER 
OF POINTS — RANGE OF EFFECTIVE NUMBER NORMALIZED 
FREQUENCY = SMOOTHING BANDWIDTH OF DEGREES STANDARD 
CASE AVERAGED (HZ) (HZ) OF FREEDOM ERROR 
4-8192 1 0-500. 0 122 8 500 
4-8192 1 0-5.98 122 8 500 
10 5.98-500.0 1.220 80 158 
4-8192 1 0-1.09 i bw! 8 500 
3 1.09-5.98 . 366 24 289 
10 5.98-500.0 1.220 80 .158 


Some additional guidelines which may be used in the selection of 


“how many “segment s>to average may be obtained from SLE EES NI 
—rreorrree 


the autocorrelation function obtained from the _segment ee 
In order to compute an inverse fourier transform, the smoothed spectra 
must consist of equally spaced frequency increments. Generally the 


spectra to be used will only be segment averaged and not frequency 
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averaged in order to preserve equal frequency spacing. In all of the 
cases which will be discussed, the segment averaged spectra was trans- 
formed using the techniques outlined in section 3. Figure 8 is a 

plot of the autocorrelation functions obtained from the spectra shown 
in Figure 6. In all cases the plots are quite similar up to a lag time 
of .2 seconds. For longer lag times there is more difference evident. 


AS the number of segments used in the frequency averaging increases,. 


the value of the autocorrelation stays closer to zero for lag times 


¢ 


from .2 to 1.0 seconds. Table 5 shows the areas of the autocorrelation 
function up to the first zero crossing and also from 0 to 4.096 seconds. 
There is up to a 25 percent difference in the area to the first zerg 
_ crossing between the different cases although the spectra of Figure 6 


appear to be virtually identical. The areas computed over the full 


range of the autocorrelation are at least one order of magnitude less 
“ ————— —~ $$ ae eisai | ——— — s 


a 


than the areas to the first zero crossing. A more detailed discussion 
of the reason for the difference in the areas is presented in the fol- 
lowing paragraphs. These two problems represent a significant mee 


if one is interested in computing an integral scale. 
: s 5 ied ties, ese 


TABLE 5 - EFFECT OF SEGMENT AVERAGING ON THE AREA UNDER 
THE AUTOCORRELATION FUNCTION 


AREA AREA 

TO FIRST 0-4.096 

CASE ZERO CROSSING SECONDS 
4-8192 0405 00142 
8-8192 0863 00157 
16-8192 0336 00195 


64-8192 .0280 .00187 
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In order to try to get a more accurate calculation of the 
autocorrelation function, a program was written to calculate the auto- 
correlation directly from the data record. This is a much more 
expensive method than the FFT technique and not as many cases were 
run. A comparison of the autocorrelations obtained using a direct 
calculation and using an inverse fourier transform of a spectra is 
shown in Figure 9. The plots with the title PROGRAM ACR were computed 
directly using a data record of the indicated length in 8 second seg- 
ments. For a 32 second record, 4 separate autocorrelations were com- 
puted and averaged in a manner analogous to segment averaging of the 
spectra. The cost of calculation was such that in the direct case, the 
computation was only carried out to a lag time of .9 seconds. Therefore, 
the only direct comparison which can be made between the plots is the 
area up to the first zero crossing. For the (32 second record, the 
area to the first zero crossing is .0353 for the direct calculation 
and .0405 for the FFT calculation. For the 64 second record, the area 
is .0362 for the direct calculation and .0363 for the FFT calculation. 
It is interesting to note the comparison in cost to obtain an auto- 
correlation via the direct method with that for the FFT technique. For 
the lower two plots of Figure 10, both of which represent a data record 
of approximately 64 seconds of real time, the direct calculation for 


100 values of time lag cost $28.00 while the FFT calculation costs 


$5.40 for 4096 values of time lag. This is a factor of 5 differences 


a 


in cost for 40 times fewer correlation points. The FFT technique also 


provides a spectrum for the cost indicated. 


(in the course of the direct calculation of the autocorrelation 


function, an interesting effect of the length of the record used in the 
———_— Eee : ; 


————————— 
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calculation was observed. The direct calculation was initially carried 
out using a record length of 2000 (2 seconds) and the maximum lag com- 
puted eomuensentad to 1000 data values (1 second). Two examples of this 
calculation are shown in Figure 10. In both cases where records of 
2 seconds each were used, the autocorrelation is negative from a time 

_ Tag of 0.3 Seepnes UO a. Cass ise oF 2.0 secands. This was not the case 
in any of the computations which used the FFT. In order to see what 
effect record length had on this negative region, the direct calcula- 
tion program was modified to use a record length of 8000 (8 seconds). 
The results of these computations for the same total length of data are 
also shown in Figure 10. The negative region from .3 to 1.0 is no 
longer predominant, and these results agree well with the autocorrela- 
tions obtained from the FFT routines as shown in Figure 9. 

An_explanation for this difference can be made based on physical 

es A time lag of .5 seconds corresponds to a frequency of 


2HZ. In a 2 second record there would only be 4 cycles at this fre- 


quency and fewer cycles at any lower frequency (longer time lags). It 


seems that 4 cycles are not enough to adequately average in the calcula- 
tion of an autocorrelation. By using a record length of 8 seconds, 
there will be 16 cycles of a 2HZ signal in one record, and the resolu- 
tion at lag times of .5 seconds will be better. Based on a limited 
_least 8 cycles of a particular frequency should be present to obtain 
adequate resolution in an autocorrelation function at_a lag time 
An additional effect of interest also arose in one case. A 


digital data tape was used which had more than one channel of data. 
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For a small portion of one record, the channels were reversed and the 
effect on the power spectral density is shown in Figure 1l. The noise 
a 

in the high frequency portion of the spectra is due to the channel 
switch. The second plot is of the same data but avoiding the record with 
the channel switch. 

The cost of the various cases run with program SEGEMNT are listed 
in Table 6. These include the computation of a power spectral density, 


an autocorrelation and plots of both using the U200 plotting routines 


available at the Engineering Research Center, Colorado State University. 


TABLE 6 - COST FOR TEST CASES - PROGRAM SEGEMNT 
CENTRAL PROCESSOR COST = $290/hr 


TIME OF 

NUMBER OF TOTAL AMOUNT 

SEGMENTS OF DATA COST 
OF LENGTH 8192 (SECONDS ) $ 

4 o2.. 107 4.00 
8 65.54 9.40 
16 131.07 $...92 
64 524..29 271 ibe 


It is important to bear in mind that all of the examples in this 
section have been calculated using a record of pressure data obtained 
using a linear transducer. Non-linear transducers or signals of a 
different type which require different frequency range or which were 
taken at a different sample rate would alter the cost figures. As such, 
these examples should only be considered as guidelines in selecting 


a scheme for digital analysis. 
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4.2 Calculation of Power Spectral Densities Using an External Core , 

FFT Algorithm 

In some applications, it is desirable to have a greater frequency 
range of the power spectral density, or resolution of the autocorrelation 
at relatively large lag times. In order to obtain either of these 
results, a long record of data must be used for each segment. In 
order to stay within the present available core of the CDC 6400, 
(140000,), the longest data record which is a power of two which may 
be used is 8192. A technique is available which allows longer data 
records to be considered by making use of disc storage and performing 
the FFT in pieces. The details of the algorithm are described by 
Brenner (3). A program titled FOR2D is available from the IBM. 
Contributed Program Library (#360D-13.4006) . This program was written 
by Norman Brenner and uses the algorithm of reference 3. The program 
allows record lengths limited only by the disc storage available on the 
computer system in use (presently between 2,000,000 and 3,000,000 for 
the CSU CDC 6400 system). This capability allows very long record 
lengths to be used if necessary. The cost of the calculations becomes 
large as longer records are used and in many cases becomes a limiting 

In order to use an external core type of program, the input data 
record is broken into a series of equal length records. It 1s necessary 
to be able to store 3 of these records in the core of the computer at 
any given time. This requirement will set the length of this array. 
The input data record is then stored on the disc and the FFT routine 


only calls a portion of the record at a time. It is important to 
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understand that this is not a segment averaged technique, but that 


the resulting sequence of points is the same that would be obtained if 


the entire data record were transformed using a computer with a very 
large core. 

A listing of a program written to utilize the external core 
technique, EXTCORE is in Appendix B3. A listing of subroutine FOR2D 
is in Appendix A2. 

The steps necessary to calculate a power spectral density are 
basically the same as were listed in section 4.1. The only differences 
between program EXTCORE and SEGEMNT are in the input and averaging. 
These differences will be pointed out with reference to line numbers 
in Appendix B3. 

In lines 145-170, the data is read from the data tape (tape 1) in 
units compatible with the length of the records to be stored in mass 
storage. These records are available to the program by calling sub- 
routine DREAD. Lines 187-205 remove the mean from the data and taper 
the data. FOR2D is called in line 209. The remainder of the program 
involves frequency averaging, output, and plotting. 

The frequency averaging is similar to that described in the 
previous section except that even the low frequency portion of the 
spectrum is frequency averaged. Since only one segment is run, some 
frequency averaging is necessary even in the low frequency portions 
in order to obtain acceptable levels of statistical reliability. 

The power spectral densities obtained from four different cases 
using program EXTCORE are shown in Figure 12. The notation in the 
figures indicates how many portions were used to make up the entire 


record. The figure in the bottom right utilized a record made up of 


ZA 


512 parts, each consisting of 1024 data elements. The averaging in 

the three shorter cases was such that the bandwidths for all three 

were the same. The fourth case (512-1024) used a different scheme of 
frequency averaging. The details of the frequency averaging along with 
the number of degrees of freedom, and the normalized standard error are 
shown in Table 7. This table can be compared with Tables 2 and 4 of 
section 4.1. It can be easily seen that as the normalized standard 


error decreases, the power spectral density function becomes smoother. 


TABLE 7 - EFFECT OF RECORD LENGTH ON POWER SPECTRAL DENSITY 
ESTIMATES , EXTERNAL CORE FFT 


NUMBER 
OF POINTS RANGE OF EFFECTIVE NUMBER OF NORMALIZED 
FREQUENCY SMOOTHING BANDWIDTH DEGREES OF STANDARD 


CASE AVERAGED (HZ) (HZ) FREEDOM ERROR 
32-1024 8 0-31.25 244 16 353 
(32.77 SEC) 16 31.25-62.50 .488 32 ,250 

128 62.50-500.0 3.906 256 088 
64-1024 16 0-15.63 .244 32 250 
(65.54 SEC) 32 15.63-31.25  .488 64 APF 
256 31.25-500.0 3.906 512 063 
128-1024 32 <7. 244 64 Te 
(131.07 SEC) 64 7.81-15.63  .488 128 125 
512 15.63-500.0 3.906 1024 044 
512-1024 16 0-1.95 030 32 250 
64 1.95-21.48 122 128 © .125 


(524.29 SEC) 
256 21.48-41.02 .488 512 .063 


a 41.02-500.0 <OTd 1024 .044 


Pi 


Table 8 shows the values of the areas under the spectra of 
Figure 12. This table is comparable to Table 3 of section 4.1. The 
first case (32-1024) shows more variation than any of the other cases 
in Table 3 or Table 8, but for many applications this error would be 


acceptable. 


TABLE 8 - COMPARISON OF INTEGRATED PROPERTIES OF POWER SPECTRAL DENSITY 
ESTIMATES, EXTERNAL CORE FFT 


CASE F(n)dn 

O 
32-1024 .978 
64-1024 1.016 
128-124 i..003 
512-1024 <997 


Correlation functions were computed for three of the example cases 
and are shown in Figure 13 along with one case calculated with program 
SEGEMNT. A trend can be seen in these figures which is similar to that 
of Figure 8. As the length of record increases, the correlation at 
larger lag times is more nearly zero. The correlations computed using 
program EXTCORE all have very much larger record lengths than those 
computed using program SEGEMNT. It would be expected that the EXTCORE 
correlations would be valid for longer lag times. A listing of the 
areas for the three correlations computed is shown in Table’9. In all 
of the cases, the area to the first zero crossing is comparable to that 
for the entire autocorrelation (0-4.096 sec). This agreement is in 


contrast with the cases shown in Table S for the shorter records of 
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program SEGEMNT. This is another example of the effect of record 

length on the calculation of autocorrelation functions. There is 

fair agreement between the areas out to the first zero crossing in 

both cases, and this may be an appropriate choice of area when only a 
limited record length is available. Care must be used in using the 

area to the first zero crossing, since not all correlations remain 

as close to zero as this demonstration case for regions beyond the first 


zero crossing. 


TABLE 9 - EFFECT OF RECORD LENGTH ON THE AREA UNDER THE 
AUTOCORRELATION FUNCTION, EXTERNAL CORE FFT 


RECORD AREA TO 
LENGTH FIRST ZERO AREA 
CASE SECONDS CROSSING 0-4.096 SEC 
32-1024 he oe 0362 . 0339 
64-1024 65.54 0378 35.76 
128-1024 Lo)<.07 . 0333 .0350 


The costs for the EXTCORE examples are listed in Table 10. These 


include the cost of all calculations and plotting. 


TABLE 10 - COST FOR TEST CASES PROGRAM EXTCORE. 
CENTRAL PROCESSOR COST = $290/hr 


DATA LENGTH COST 

CASE SECONDS $ 
32-1024 32:77 12.10 
64-1024 65.54 22.59 
128-1024 131.07 46.37 
512-1024 524.29 210.30 
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4.3 Comparison of Program SEGEMNT and Program EXTCORE 

Many of the differences between and advantages of segment averaging 
and external core approaches are apparent after reading the previous 
section. These differences and advantages will be briefly summarized 
in order to point out the most significant. 

The major advantages of the external core technique are that it 
allows a greater frequency range in the spectrum and provides an 
autocorrelation function which is valid at relatively long lag times. 
The advantage of being able to obtain more points at low frequency in 
the spectrum is offset somewhat by the need to perform some type of 
smoothing in order to obtain a statistically reliable value. For the 
external core case, the smoothing will be accomplished using frequency 
averaging which will reduce the number of data points available at the 
low frequency end of the spectrum. 

The autocorrelation which may be obtained using the external 
core technique is of higher quality at higher lag times than the 
autocorrelation which may be obtained using segment averaging. This 
increase in quality is obtained at a corresponding increase in cost 
of computer time. This extra cost may be necessary if an accurate 
measure of integral scale is desired. The integral time scale and the 
low frequency end of the power spectral density are directly related, 


co 


(F(0) = 4 f R(t)dt), and if the low frequency end of the spectra has 
O 

a standard error of .5, there can be up to 50 percent error in the 

integral scale. 


The major advantage of segment averaging is cost. In all cases, 


comparable quality power spectral densities can be obtained (based on 
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normalized standard error) for from 1/3 to 1/8 the cost using segment 
averaging instead of external core techniques. 

The core requirements for each case are comparable based on the 
array sizes used in the example programs. Changing array sizes in 
either of the programs would have an effect on the core required, but a 
comparable change would have to be made to both programs and the core 
requirements would still be comparable. 

A general guideline in selecting a technique would be to use 
segment averaging unless a special requirement exists which requires 


the external core technique. 


D TWO CHANNEL CALCULATIONS--CROSS-SPECTRAL DENSITIES AND CROSS- 
CORRELATIONS 


Some applications require information concerning the relationship 
between two time series in either the frequency or time domains. Once 
the techniques described in the previous sections are understood, the 
computation of functions describing these relationships can be readily 
accomplished. Most computations of multichannel functions begin with 
a cross-spectral density function, a complex quantity. Once the cross- 
spectral density function is obtained, a number of additional quantities 
can be computed. A brief discussion of some of these functions can be 
found in Bendat and Piersol (1), pp. 25-34. 

The equation for the cross-spectral density of two time series 
x(t) and y(t) is given by the equation G y(t) = = x*(£.) Y (f)- 
(x* is the complex conjugate of the transform of the x(t) time series.) 
Thus, once the transforms of two simultaneous time series are available, 
the cross spectral density, and any other related quantities may be 


computed. As brief examples of both computation and averaging 
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techniques, programs which compute a coherence function and a 
cross-correlation coefficient will be discussed. 

The two most important aspects of these programs are the 
techniques of data storage and averaging. The data storage is common 
to both programs and will be explained with reference to PROGRAM 
CSPECT2 (Appendix B4). The single channel transforms have been com- 
puted and are stored on a master data tape (tape 2) as separate files, 
with each segment a separate record (logical record) of the file. The 
input portion of the program (lines 90-105) reads each file from 
tape 2 and stores them on tape 3 and tape 4 for X(f,) and Y(f) 
respectively. All reads and writes are done using unformatted binary 
reads and writes. The use of this type statement instead of a format- 
ted read or write results in savings of from 50 percent to 90 percent 
in the required computer central processor costs. 

As shown in previous sections, some method of averaging will be 
required to obtain statistically reliable estimates. In the examples 
segment averaging is used as the primary method. It is necessary to 
segment average the cross-spectral density function and not the single 
channel transforms. Therefore, in lines 110 and 111 the single channel 
transform for each segment is read and the segment averaged cross- 
spectral density function is computed. In the calculation of the 
coherence function (lines 113-130) the cross-spectral density is also 
frequency averaged prior to the final calculation of the coherence 
function (lines 138-156). 

A second example program which calculates a cross-correlation 
function (PROGRAM CSPECT3) is listed in Appendix BS. The input and 


smoothing sections of this program are the same as those in CSPECT2. 
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The cross-correlation is obtained from an inverse fourier transform of 
the cross-spectral density function. A segment averaged cross-spectral 
density 1s computed in lines 106-113 and reflected in lines 117-122. 
The cross-spectral density is reflected such that the real part is 

an even function and the imaginary part an odd function. The reflected 
cross-spectral density is transformed in line 126 to obtain a cross- 
correlation coefficient. The cross-correlation coefficient can be 
calculated directly from the time series, and a comparison of a direct 
computation and a FFT computation is shown in Figure 14. The two 
results are virtually identical. 

This brief section shows just two of the many cross-channel 
computations possible. Costs of the different calculations will vary 
with the application and no definite guidelines can be stated. The 
two most important aspects of cross-channel calculations are 
(1) use of binary write and read statements (2) averaging the cross- 


spectrum and not the single channel transforms. 
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